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ABSTRACT 

Given the present distribution of mass tracing objects in an expanding universe, we 
develop and test a fast method for recovering their past orbits using the least action 
principle. In this method, termed FAM for Fast Action Minimization, the orbits are 
expanded in a set of orthogonal time-base functions satisfying the appropriate bound- 
ary conditions at the initial and final times. The conjugate gradient method is applied 
to locate the extremum of the action in the space of the expansion coefficients of 
the orbits. The TREECODE gravity solver routine is used for computing the gravi- 
tational fields appearing in the action and its gradients. The time integration of the 
Lagrangian is done using Gaussian quadratures. FAM allows us to increase the number 
of galaxies used in previous numerical action principle implementations by more than 
one order of magnitude. For example, orbits for the ^ 15, 000 IRAS PSCz galaxies 
can be recovered in 12, 000 CPU seconds on a 400MHz DEC- Alpha machine. FAM 
can recover the present peculiar velocities of particles and the initial fluctuations field. 
It successfully recovers the flow field down to clusters scales, where deviations of the 
flow from the Zel'dovich solution are significant. We also show how to recover orbits 
from the present distribution of objects as in redshift space by direct minimization of 
a modified action, without iterating the solution between real and redshift spaces. 

Key words: cosmology: theory - gravitation - dark matter - large scale structure 
of Universe 



1 INTRODUCTION 



In the standard cosmological paradigm, the present distribution of galaxies and their peculiar motions are the result of 
gravitational amplification of tiny initial density fluctuations. Accordingly, the mildly non-linear large scale structure observed 
today contains valuable information on the initial fluctuations. Given either the present galaxy distribution or the peculiar 
velocity field one can recover the growing mode of the initial density field (Peebles 1989, Weinberg 1992, Nusser & Dekel 
1992, Gramman 1993, Giavalisco et al 1993, Croft & Gaztanaga 1998, Narayanan & Weinberg 1998). Methods for recovering 
the initial fiuctuations can be very rewarding as one can directly address statistical properties of these fiuctuations (Nusser & 
Dekel 1993). As has been shown by Nusser & Dekel (1992,1993) the initial density field recovered from peculiar velocity fields 
has a strong dependence on the value of the density parameter, fl. On the other hand, a recovery from the galaxy distribution 
depends very weakly on SI. Therefore by matching the statistical properties from the two recovered initial density fields one 
could provide an estimate for fl. Gravitational instability theory, also provides tight relations between the present density and 
peculiar velocity fields (e.g., Nusser et al 1991). These relations have been an important tool in large scale structure studies, 
in particular for model independent estimates of the cosmological parameters. For example, a comparison of the measured 
peculiar velocities with those predicted from galaxy redshift surveys can yield Q. (Davis, Nusser & Willick 1996, Willick et 
al 1996, da Costa et al 1998, Sigad et al 1998, Branchini et al 1999). These comparisons are done under an assumed form for 
the biasing relation between the mass and galaxy distribution. So deviations from the theoretical velocity-density relations 
may serve as an indication to the way galaxies trace the mass and hence to the interplay between galaxy formation and the 
large scale environment. 
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Most of the promising metliods for recovering the initial growing mode and for relating the present peculiar velocity and 
density fields, are based on the least action principle (LAP). Peebles (1989) (hereafter P89) has pointed out that the equations 
of motion can be derived from the stationary variations of the action with respect to orbits subject to fixed final positions and 
vanishing initial peculiar velocities. P89 proposed minimizing the action with respect to the coefficients of an expansion of 
the orbits in terms of time-dependent functions satisfying the appropriate boundary conditions. Methods based on LAP are 
very powerful as the true orbits arc recovered to any accuracy depending on the number of functions used in the expansion, 
at least in the laminar flow regime. They also provide simultaneous estimates of the present peculiar velocities and the initial 
fluctuations from a given distribution of galaxies. 

So far, the action principle has been applied to study the dynamics of the Local Group of galaxies (P89, Peebles 1991, 

1994, Dunn & Laflamme 1993, Branchini & Carlberg 1994) and to recover the peculiar velocity from the distribution of about 
1100 galax;ies within a redshift of 3000 km/s (Shaya, Peebles & Tully 1995). The application of LAP to larger data sets have 
been hindered by the heavy computational burden needed in current methods for minimizing the action. The LAP is the only 
technique with which one can probe the nonlinear behavior to any accuracy as opposed to Zel'dovich type approximations 
which break down near large mass concentrations like clusters of galaxies. Therefore, a fast method for implementing LAP 
can be rewarding. Such a method would enable us to apply LAP to large data sets like the various IRAS galaxy redshift 
surveys (1.2 ,]y. Fisher et al 1995a, and PSCz, Saunders 1996) and portions of the Sloan Digital Sky Survey (Gunn & Knapp 
1993). In contrast to previous numerical implementations of LAP, which use direct summation over inter-particle forces, 
we propose here an efficient method based on the TREECODE to compute gravitational forces and the potential energy. 
Also, our method expands the orbits in orthogonal time base functions and perform the time integration using the Gaussian 
quadratures method. All this increases the speed of the calculation by more than one order of magnitude over any previous 
implementation of LAP. We refer to our method by the acronym FAM for Fast Action Minimization. 

The outline of the paper is as follows. In section 2 we review the least action principle and describe FAM. In section 3 we 

demonstrate the robustness of and show tests of the recovered peculiar velocities. In section 4 we describe extensions of the 
FAM to flux limited surveys, application from redshift surveys and biased distribution of galax;ies. We conclude in section 5 
with a discussion of our results and possible applications of FAM. 



2 THE BOUNDARY VALUE PROBLEM AND THE LEAST ACTION PRINCIPLE 

We follow the standard notation in which a{t) is the scale factor, H{t) = a/a is the time-dependent Hubble factor, Q = pb/pc 
is the ratio of the background density, pt, of the universe to the critical density, pc = 3H^ /8ivG. 

Assume that the underlying mass density field in a spherical volume V is sampled, in an unbiased way, by a discrete 
distribution of N gala:xies (particles). In this sampling, if the average mass density in any cell of volume (5K « V is p^ 
then the number of particles in that cell is drawn from a Poisson distribution with mean n{p^ / pi,)SV , where n = N/V is 
the mean number density over the large volume V and we have assumed that the average mass density in V is pt- Instead of 
using the usual time variable, t, we describe the evolution of the system in terms of the linear growing mode, D{t), of density 
perturbations (e.g., Peebles 1980). Let Xi denote the comoving coordinate of the i**^ particle and 0i = dXi/dD its velocity 
with respect to the time variable D. Neglecting interactions between matter interior and exterior to the volume V, system of 
particles obeys the the following Euler equations 

dOi 3 13 1 n 

dD^ 2d'^'~ 2D^ Pin)^^^'^' ^ ' 

where /(f2) = dlnD/dlno « OP'^ is the linear growth factor (e.g., Peebles 1980) and g represents the peculiar gravitational 
force field per unit mass. These equations are supplemented by the Poisson equation 

V ■ gix) = -Six) , (2) 

which relates the divergence of g to the mass density contrast 5{x) = p{x)/pb — 1 at any point X in comoving coordinate space. 
We can approximate 6 from the discrete particle distribution by 



S{x) = yJ2s''{x-Xi)-l, (3) 
where 5° is the Dirac delta function with unit integral over the volume V. Therefore, the field g is given by 
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(4) 



Equations (^) and ^ constitute the equations of motion governing the evolution of the system of particles. We do not include 
the continuity equation in the equations of motion. The continuity equation is a constraint equation obeyed automatically 
as the particles move according to and (^. The equations of motion involve second order time derivatives. Solving these 
equations can be seen either as an initial value problem or as a boundary value problem (Giavalisco et. al. 1993). Numerical 
N-body codes (e.g., Hockney & Eastwood 1981) are the usual tool for solving the relevant initial value problem where the 
positions and velocities of particles are specified at a given time. But, recovering the orbits of particles given their present 
positions is a boundary value problem where the solution must yield a homogeneous distribution of particles at the initial time 
D — 0. P89 stated the cosmological boundary value problem in the context of the least action principle. P89 also suggested an 
approximation to the orbits by minimization of the action with respect to a particular choice trial functions. In our notation, 
the action of the system of particles is 



S = 



2 3 1 

' + 



2 Di/2 f2(n) 



Aim ^-^ \Xi — Xj 



+ 



(5) 



where we have arbitrarily set D = 1 at the present time. Following P89 we minimize the action with respect to the coefficients 
of an expansion of the orbits by means of time dependent base functions {qn{D),n = 1 • ■ ■ Umax}- We write the position of 
each particle Xi{D) for D < 1 as 



Xi{D) = XiM + ^ qn{D)Ci, 



(6) 



where Xi,o is the position of the particle at D = 1 and the vectors Ci,n are the expansion coefficients with respect to which 
the action is to be minimized. Since Xi{D = 1) = Xi,o, we choose the functions qn{D) such that qn{D = 1) = 0. By taking 
derivatives of (pi) with respect to D, we find that the velocity, di, is 



^^{D) = p„{D)Ci,n, 



(7) 



where p„ — dq„/dD. Stationary variations of the action with respect to d yield the following set of equations 
dS 



/" J 7-1 7-i3/2 fl , 3 ^ 



0. 



(8) 



If we integrate by parts the term involving the velocity in the previous equation, we arrive at 



lim {D''''^q„d, 



»1 



dBj 3_1_ _ 3 1 n 

dD^ 2D " 2D2 p(Q.)^' 



0. 



(9) 



Without the boundary terms on the left, these equations are equivalent to the equations of motion averaged over time with 
weight functions D^^^qn. The boundary terms are individually eliminated by the imposing the following two constraints on 
the time-base functions (Peebles 1989), 



g„(D = 1) =0 



and 



D^''\n(D)e(_D) = 0. 



(10) 



2.1 The homogeneity condition and the time-base functions 

We will work directly with the functions p„ = dqn/dD rather than g„. The first constraint in ( |lo|) means that the positions 
of particles at _D = 1 remain unchanged when varying d^n- The second constraint is very fiexible, it merely implies that 
dlnp„(D)/dln_D > —5/4 for D ^ 1. However the functions p„ must lead to homogeneous initial particle distribution. A 
sufficient condition for the homogeneity is that the time dependence of the particle velocities near D = matches that of 
the linear velocity growing mode which, with respect to the time D, is a constant. Orbits with velocities with initial time 
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dependence like that of the decaying mode (cx _D~^''^), do not necessarily lead to homogeneous initial distributions, although 
the decaying mode is derived under the assumption of small perturbations. Therefore, initial homogeneity implies that one 
of the Pn must be a constant and the rest increasing functions of D. There are many functions satisfying our boundary 
conditions. Here we choose pn to be linear combinations of 1, D, ■ ■ ■ _D"™°^ (Giavalisco et. al. 1993) which satisfy the 
following orthonormality condition 



[ dDD''^^p„{D)p^{D) = 5^,„ 
Jo 



(11) 



where 5^ is the Kronecker delta function. Orthonormality will prove useful in the numerical minimization of the action using 
the conjugate gradient method. The functions p„ can be constructed using the Gramm-Schmidt algorithm, however, in this 
case they can be derived from the expression 

P^iD) = A„-^— [D'/'D-il - DT] , (12) 

where An are normalizing constant. With the orthonormality condition ( pT| ) the expression for the action gradients dS/dC" 
(|) becomes 

=C.,n + U'l^?,9.iD)dD. (13) 



dCi,„ 2 J.. p 



2.2 The numerical action minimization 

Given the orbit expansion ^ and the base-functions (|l^, the problem of recovering the orbits reduces to finding the value 
of the coefficients Ci,„ where the action ^ has a minimum. Our method for minimizing the action, FAM, is based on the 
conjugate gradient method (CGM) (e.g., Press et. al. 1992). An efficient implementation of CGM calls for a fast way of 
computing the action and its gradients with respect to C. Most of the computational cost comes from the potential energy 
term in the action and the gravitational force field, g, in the gradients (|l3|). These quantities involve summation over pairs 
and they have to be computed various times, in each step taken by CGM, for an accurate estimate of their integrals. We can 
achieve a significant improvement over previous schemes for minimizing the action if instead of computing the gravitational 
forces by direct summation, we use the TREECODE technique (e.g., Bouchet & Hernquist 1988). Although any of the 
fast techniques like the particle-mesh (PM) or particle-particle particle-mesh (P'^M) (Hockney & Eastwood 1981) or the 
adaptive P^M (Couchman 1991) can be used, the TREECODE is particularly suitable for our purposes since it can readily 
be implemented for particle distributions in a spherical region as in whole sky galaxy catalogs. We reduce even further the 
required number of gravitational field calculation by performing the time integration in the expression for the action and its 
gradients using the Gaussian quadrature scheme with D^^^^ weights (Giavalisco et al 1993). 

The boundary value problem is almost certain to have more than one solution (P89, Giavalisco et. al. 1993). Part of 
the reason for this is that one of the boundary conditions only prescribes time dependence of the velocities near the initial 
time D = and does not specify their amplitude. This is not sufficient for a unique solution, in the presence of orbit mixing 
regions. So the action can have many minima corresponding to different solutions of the time averaged equations of motion. 
Therefore, we expect the minimum found by CGM to depend on the initial guess. However, our purpose is to recover motions 
on large scales which means that out of all minima we would like to find the one corresponding to orbits which do not deviate 
significantly from the Zel'dovich straight line approximation. A reasonable choice for the initial guess is then: d^n ~ for 
n > 1, and Ci,i obtained from ( |l3| ) by substituting gi{D) — Dg^ q where g^ Q is the gravitational force field at D = 1. This 
means that the initial guess for the orbit of each particle is a motion in straight line with velocity given by the gravitational 
force field according to linear theory. 



2.3 The dependence on Ho and Q 

We would like to clarify here the dependence of the equations of motion and the recovered orbits on the present value of 
Hubble factor. Ho, and the density parameter, il. Since _D is a dimensionless variable, the spatial coordinate x and the velocity 
d = Ax/ AD have the same units, e.g., Mpc. We could also work in km/s if we define HqX to be the spatial coordinate. Neither 
the action nor its gradients (^) involve the Hubble factor, so the recovered orbits are completely independent of the units 
with which we choose to express the orbits. 
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Figure 1. Recovered orbits from tlie FAM. Filled dots show present time positions for a random selection of N-body particles contained 
within a slice of thickens 10 Mpc. The solid lines represent their projected orbits. The upper plot is an enlargement of the central 
region in the lower plot. 

The equations of motion expressed in terms of the time variable, D, are almost independent of Q and the cosmological 
constant (e.g., Weinberg & Gunn 1990, Nusser & Colberg 1997, Mancinelli & Yahil 1995, Gramman 1993). This is because 
of the weak dependence of f2//^ ~ Jl""'^ on the cosmological parameters in 0. However, as we shall see in subsection 4.2, the 
orbits recovered from the distribution of galaxies in redshift space rather in real space, will have a non-negligible dependence 
on SI. One could then obtain orbits in a f2 7^ 1 universe by appropriate scaling those recovered assuming flat universe. 



3 TESTING FAM WITH AN N-BODY SIMULATION 

FAM involves a number of parameters. These include the force softening scale in the TREECODE gravity solver, the number 
of time-base functions, the convergence tolerance parameter in CGM, (see Press et. al. 1992 for details). In addition, we have 
the initial guess Ci,n required by CGM. We resort to an N-body simulation to demonstrate the robustness of the method 
to these parameters and the initial guess. We use a high resolution simulation (Cole et. al. 1998) of cold dark matter in a 
flat universe with a cosmological constant. The matter density parameter at the final output of the simulation is Qo = 0.3 
The simulation contained 192'^ particles in a periodic cube of side 345.6 Mpc. The simulation is normalized such that 
at the final output the linearly extrapolated rms value of the density fluctuations in a sphere of Sh^^Mpc is as = 1.13. 
The simulation was run using a modified version of Couchman (1991) AP^M N-body code with a force softening parameter 
of 0.27 h-^ Mpc. We test FAM using the distribution of 1.5(T0*) particles selected at random from a spherical region of 
radius 80 Mpc in the simulation. We will refer to a standard FAM recovery as the one in which the orbits are expanded 
in six time-base functions given by (|l^, the tolerance parameter in CGM is 10"'*, and the initial guess for d.n is given 
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Figure 2. Robustness tests of FAM. To the left: each panel shows velocities (one component) ' recovered with Ci^n = as the initial 
guess vs those recovered with standard FAM. The panels correspond, respectively, to different values for the force softening parameter. 
To the right: each panel compares velocities recovered with two values of the convergence tolerance parameter as indicated at the top. 



from linear theory, as described at the end of subsection (2.2). The total number of free coefficients Ci,n in standard FAM is 
3 X 1.5(T0'') X 6 = 2.7(-10^), including the 3 spatial components. 

We will focus on the performance of FAM at recovering the particle velocities at the final output of the simulation. 
However, it is instructive to visually examine the recovered orbits. The solid lines in figure 1 are two dimensional streamlines 
of particles contained, at the final time, in a slice of thickness of 10/i~^ Mpc. The dots indicate the present positions of the 
particles. The streamlines shown in the plot correspond to orbits recovered with standard FAM and a force softening parameter 
of 0.5/i~^Mpc. The deviations of the streamlines from straight lines are significant, especially in high density regions. These 
deviations are an indication of the failure of the Zel'dovich approximation. 

We first assess the robustness of the recovered velocities against changes in the initial guess for d^n- To do that we have 
compared the solution of standard FAM with the solution obtained = as the initial guess. The three panels 

on the left hand side in figure 2 compare between the FAM recovered velocities in the two casesQ. The three panels show, 
respectively, results for three values of the force softening parameter, as indicated in each panel. Changing the initial guess did 
not introduce any systematic differences in the recovered velocities for all three values of the softening parameter. The scatter, 
although not negligible, is small compared to the scatter between the recovered true velocities (see figure 4 below). The right 
hand side of figure 2 compares the standard FAM velocities with those recovered with a convergence tolerance parameter of 
10~^. The correlation between the two velocities is very tight and no systematic differences are detected. 

Having established the robustness of FAM we now proceed to check how well it reproduces true velocities of particles 
in the simulation. In figure 3 we show the scatter plots of recovered vs true velocities of randomly chosen particles. In 



* In this and the following figures we plot comoving peculiar velocities V = da;/dt = Hof{Qo)0 measured in units of km/s 
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Figure 3. Recovered vs true velocities for random particles in the simulation. Left, middle and right columns correspond, respectively, 
to recovery from, standard FAM, Zel'dovich approximation, and linear theory. Shown are velocities smoothed with a top-hat window of 
radius 5 (top row), 3 (middle), and 1 Mpc (bottom). Displayed in each panel, the parameters of the linear regression of recovered 
on true velocities. The 45° solid lines are drawn to guide the eye. 

addition to velocities recovered with standard FAM (left column) we show results from the Zel'dovich approximation (middle 
column) and linear theory (right column). Velocities in the Zel'dovich approximation were obtained by running FAM with 
rimax = 1, i-e, with straight line orbits of the form Xi{D) = Xo,i + Dd^i. The Zel'dovich velocities should coincide with 
those which would have been recovered by the PIZA method of Croft & Gaztanaga (1998). The linear theory velocities are 
simply _ffo/(no)</o,i where Qq ^ is the gravitational force field obtained from the particle distribution at the final time. In all 
cases the softening parameter is 0.5/i~^Mpc in the computation of the gravity field. The recovered and true velocities in each 
panel with a top-hat window of the same length. The three rows show results for three smoothing lengths, as indicated in the 
figure. The top-hat smoothing replaces the velocity of each particle by the mean velocity of its neighbors within a distance 
equal to the smoothing width. The parameters of the best linear fit are shown on the top left corner of each plot. Linear 
theory performs very poorly, even when smoothing on scales as large as 5 . Evidently, the Zel'dovich approximation is a 
significant improvement over linear theory. Yet, a close visual inspection reveals a systematic bias in the Zel'dovich velocities. 
This is confirmed quantitatively by the parameters of the linear regression. The FAM velocities are almost unbiased even for 
the smallest smoothing width. Note also the small scatter of 150 - 250 km/s between the FAM and true velocities. 

The superiority of FAM over the Zel'dovich approximation can be further appreciated from figure 4 where we plot the 
unsmoothed velocities. For both in FAM and Zel'dovich we we show results with force softening of 0.25 Mpc (top panels) 
and 3 /i^^ Mpc (bottom panels) in the TREECODE. The lower softening matches the the one used for the force calculation in 
the original N-body simulation. The Zel'dovich solution fails to reproduce nonlinear motions for both values of the softening 
parameter. It typically underestimates the true velocities. FAM, however, provides unbiased estimates of the true velocities 
if the force softening is close to the one of the N-body simulation. FAM cannot model accurately the structure of the orbits 
on scales that the force softening length. This is the cause of the bias in FAM velocities recovered with the high value of the 
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Figure 4. FAM vs. Zel'dovich. Plotted are unsmoothed velocities. To the right: Zcl'dovich t;*. True. To the left: LAP vs. true. Shown 
are velocities recovered with softening parameters of 0.25 Mpc (top plots) and 3 Mpc (bottom). The parameters of the linear 
fit are shown in each panel 

softening parameter (bottom right panel). The random errors in the unsmoothed FAM velocities are large (~ 500km/s) but 
they go down by a factor of 2 when smoothing on a scale of 1 Mpc (see bottom left panel in figure 3). 



4 EXTENSIONS OF FAM 

So far we had in mind an ideal situation in which we have a perfect volume limited unbiased distribution of galaxies in real 
space. All sky surveys like the IRAS sample, provide us with the redshift coordinates of galaxies. They are typically flux 
limited so that the observed number density of galaxies is a decreasing function of distance. Galaxies may also be biased 
tracers of the underlying mass density field. We will now outline how FAM can deal with these issues. 



4.1 Selection effects and shot-noise 

Suppose that a flux limited catalog is obtained from a parent volume limited catalog. If the observed number density in real 
space in a fiux limited catalog is no (a;), then the corresponding number density in the volume limited catalog is no{x) / (j){\x\) , 
(j) is the selection function. The gravitational field g can then be approximated by 




(14) 
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where 0i = ^(Xi.o) and ni is an estimate of the number density in the volume hmited catalog, for example, ni ~ 
where the sum is over galaxies within a spherical region of volume V around the observer. The potential energy term in 
the expression for action is also modified accordingly. So for flux limited surveys the gravitational potential and force fields 
are computed by the TREECODE technique with each particle having a mass proportional to the inverse of the selection 
function at its final position. We now examine the error in the recovered velocities and positions of particles as a result of 
the discrete sampling of the mass density field. We define the error covariance matrix between two quantities X and Y as 
< Z\X Ai^ > where the symbol < . > denotes averaging over many different discrete samplings of the underlying mass field 
with the same selection function and number of particles as the original flux limited distribution. The A denotes the difference 
between the true and recovered values. The true value being obtained by application of LAP on a sampling with an infinite 
number of particles, but with the same selection function as in the dilute distribution. We focus on computing the velocity 
error covariance matrix. The calculation for the recovered positions is similar. Using (^), The velocity error covariance matrix 
between two particles i and j can be written in terms of the coefficients C error matrix. 



< Afli Afij ^J3„pm < AC,,„ ACj,™ > , (15) 



By setting to zero the action gradients (|13[) we find 



<AC,„AC7,„>= ^^'di3di3'^^^^ <Aff,^>, (16) 

were quantities with and without the prime symbol are evaluated at times D' and D, respectively. This expression involves the 
gravity error covariance matrix computed between errors at different times. We can estimate this matrix under the assumption 
that the deviations from the Zel'dovich solution do not affect on its value. According to the Zel'dovich approximation, the 
gravity force acting on a particle at a time D is Qi = D{f^ /Q)giQ, where g^Q is computed at the final time. With this 
approximation (^) becomes 

< ACi,„ ACj.m >= B„,m < Z^i,o ^j,0 >, (17) 

where Bn,m = (9/4) q„q'm{DD'y^^{ff')^ /{ni}')dDdD' . The calculation of the gravity error covariance at the final time is 
straightforward (Yahil et. al. 1991) and the result is 

A„ 1 1 (Xi -Xk){Xj -Xk) , . 



4.2 Redshift space distributions 

Typically galaxy surveys provide redshifts and angular positions of galaxies in the sky. Redshifts of galaxies differ from their 
distances as a result of peculiar velocities along the line of sight. This causes differences between the distribution of galaxies in 
real and redshift space (e.g., Kaiser 1987, Hamilton 1993). These differences are referred to as redshift distortions. Method for 
recovering the velocity from the redshift space distribution of galaxies in the nonlinear regime have so far relied on iterations 
between redshift and real space (Shaya et al 1995, Yahil et. al. 1991, Schmoldt and Saha, 1998). Here we show that FAM can 
be extend to treat redshift distortions by direct minimization of a modified action, without the use iterations. For simplicity, 
we restrict the analysis to volume limited surveys. We define the comoving redshift coordinate, So, of a galaxy at the present 
time, as 

So = HoXo + {Vo ■ So)so (19) 

the subscript refers to quantities at the present time and the unit vector So points in the direction of the line of sight to the 
galaxy. The comoving peculiar velocity of the galaxy is V = (dx/dt) — DH f{Q)0. 

The redshifts of galaxies are given. So the expansion of orbits in terms of time-dependent functions has to be such that 
the redshift coordinates as given by (^9|) are fixed under variations of the expansion coefficients. Defining parallel (") and 
perpendicular (^) directions to the line of sight at the present time, we write the expansion of the orbits as 

n n 
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(20) 



where po,n ~ Pn{D = 1) and Qn{D) = qn{D) — fopo,n- This expansion of the orbits satisfies the boundary conditions of fixed 
redshifts and angular positions on the sky. The role of the Hubble constant Ho is adjusting the units. The trivial dependence 
on Ho can be completely eliminated by working with HoX instead of x. Stationary first variations of the action with 
respect to subject to the boundary conditions (|l^ yield 



dOi idf 3 1 n J 
AD ^ 2"D ^ 2D2p^^ 



(21) 



which is the time averaged equation of motion of a particle in the plane perpendicular to its sight-line at the present time. 
On the other hand, stationary variations with respect to C" yield 



dDD^^^Q 



d^ 3e[_3_l__n II 
dD 2 D 2D2/2^' 



(22) 



where Qo,n = Qn{D = 1) = —foPo.n and 0q - — '^po,nCK These equations differ from the time averaged equations of motion 
by a boundary term. This term can be eliminated by adding to the action (P) a kinetic energy term corresponding to the line 
of sight parallel degree of freedom, as follows 



^-«+|E(<.)' 



(23) 



Minimization of the the modified action readily yields the orbits expansion parameters Ci,n (see Schmoldt and Saha 1998, for 
a similar treatment of this problem). 

The recovered orbits from redshift space depend on Q through the /o in the expansion ( po[ ) . The effect of this dependence 
in the linear regime is elucidated in Nusser & Davis (1994). Since linear theory overestimates the true velocities (Nusser et. 
al. 1991), this fl dependence will be weaker in the nonlinear regime. 



4.3 Biased distributions 

So far we have assumed that the number of galaxies in a small cell is proportional the the average mass density in the cell. 
However, galaxies are most likely to be biased tracers of the mass distribution as indicated by the relative bias between galaxies 
of different luminosities and morphological types (e.g, Loveday et al 1995). For simplicity of notations we discuss here how 
our scheme can incorporate biasing only for volume limited galaxy distributions in real space. Suppose we are given smooth 
versions of the galaxy number density and the mass density fields. Working with smoothed fields seems reasonable because 
galaxy formation at a given point is likely to be affected by the nearby dark matter environment. We also assume that the 
smoothing scale is large enough such that the smoothed galaxy number density field is not contaminated by shot-noise. Let 
5^ and 5 be, respectively, the galaxy number density contrast and the mass density contrast, both smoothed with the same 
smoothing window of width fixed by the physical processes involved in galaxy formation. For unbiased galaxy distribution 
5^ = S and for the familiar linear biasing 5^ = bS where b is the linear bias factor. Here 5^ is a nonlinear function of S, which 
we assume to be local and deterministic (e.g., Dekel & Lahav 1998). We characterize the biasing relation at any point in space 
by the ratio 



1 + Ss' 



Our definition of biasing in terms of smooth fields inevitably implies that the galaxy distribution does not contain information 
on the structure of mass density on scales smaller than the smoothing scale length. Given W we define the following density 
field 

g(x) = ^W(^)l V5°(x-a;.) . (24) 
n V — ' 
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This field serves as an unbiased estimate of smootlied underlying mass distribution. The gravitational force field, g, is then 
estimated from q by 



where W; = W[Xi (£*)]• If the galaxies and the mass particles share the same velocity field then the continuity equation implies 
that W remains constant along the streamlines so that >V[a5i(i3)] = W(!Bi,o) (e.g., Nusser & Davis 1994, Pry & Gaztanaga 
1995). Therefore, it is sufficient to specify the biasing relation at the present time. The net effect of biasing is that it changes the 
weight assigned to each particle in the calculation of gravitational fields. This can readily be incorporated in the TREECODE 
by assigning to each particle a mass proportional to W{xifi). 



5 DISCUSSION 

We have presented a fast method for solving the boundary value problem of recovering the orbits of particles from their 
present positions assuming homogeneous initial conditions. The method which we term FAM is based on P89 implementation 
of lest action principle in a cosmological contest. It can be applied to distribution of galaxies in redshift space. It can also very 
easily incorporate any local biasing relation. FAM is suitable for recovering orbits from large galaxy redshift surveys such as 
the PSC2:. It can also bo applied to large portions of the future Sloan Digital Sky Survey. 

We have described the method assuming that the number of expansion coefficients is the same for all galaxies. However, 
this need not be the case. For example, galaxies in low density regions can be assumed to move along straight line just like in 
the Zel'dovich approximation. This can significantly speed up FAM especially for the flux limited surveys with dilute galaxy 

distribution at largo distances from the observer. 

We have used an N-body simulation to show that FAM recovers very well the final velocities from a given volume limited 
particle distribution in real space. However, galaxy surveys provide galaxy distribution in redshift space. Redshift distortions 
introduce the nuisance of multi- valued zones where particles overlap in redshift space while they are far apart in real space. 
FAM allows a recovery of the orbits non-iteratively from redshift space data by direct minimization of a modified action. We 
believe that this should mitigate the effects of multi-values zones in the recovered orbits. Tests of FAM recovery from both 
flux and volume limited distributions in redshift space are underway. 

In this work we concentrated on how well FAM can reconstruct objects' peculiar velocities and we did not check how well 
it can recover initial fluctuations. Judged by its superiority over the Zel'dovich solution at recovering the present velocities, 
FAM is expected to perform well at recovering the density fluctuation at any time in the past. We have outlined FAM can 
incorporate possible biasing between the mass and galaxies. We have shown that we only need to specify a biasing relation at 
the present time. The recovery of the orbits is sensitive to the biasing relation. So one can tune the biasing relation so that 
the clustering amplitude of the recovered mass density field varies with time according to hierarchical structure formation. In 
hierarchical clustering, the evolution of clustering amplitude, as measured for example by the correlation function, is almost 
independent of the linear mass power spectrum (e.g., Hamilton et. al. 1991, Jain, Mo & White 1995, Peacock & Dodds 1995), 
and of the cosmological parameters if expressed as a function of the linear density growing mode, D (e.g., Nusser & Colberg 
1997). Therefore, one can determine the biasing relation independent of the cosmological parameters and the details of the 
dark matter model. 
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